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ABSTRACT 

We present an exact solution of a probabilistic cellular automaton for traffic with open 
boundary conditions, e.g. cars can enter and leave a part of a highway with certain proba- 
bilities. The model studied is the asymmetric exclusion process (ASEP) with simultaneous 
updating of all sites. It is equivalent to a special case (v max = 1) of the Nagel-Schreckenberg 
model for highway traffic, which has found many applications in real-time traffic simulations. 
The simultaneous updating induces additional strong short range correlations compared to 
other updating schemes. The stationary state is written in terms of a matrix product so- 
lution. The corresponding algebra, which expresses a system-size recursion relation for the 
weights of the configurations, is quartic, in contrast to previous cases, in which the algebra 
is quadratic. We derive the phase diagram and compute various properties such as density 
profiles, two point functions and the fluctuations in the number of particles (cars) in the 
system. The current and the density profiles can be mapped onto the ASEP with other 
time discrete updating procedures. Through use of this mapping, our results also give new 
results for these models. 
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1 Introduction 



In this paper we study a simple probabilistic cellular automaton which describes the flow — 
of particles, automobiles, or some other conserved quantity — through a one dimensional 
system. The particles (or cars) of the model move on a finite lattice; at integer times they 
simultaneously attempt to hop one site forward, succeeding with probability p if the site in 
front of them is empty. We are interested in the case of open boundary conditions (OBC), 
in which, simultaneously with the hopping of particles along the lattice, a particle enters 
the system with probability a at the leftmost site if that site is empty, and if the rightmost 
site is occupied then the particle on that site exits with probability (3. We use a matrix 
product ansatz to give a complete solution of this model. 

In recent years, cellular automata models for traffic flow have gained much attention, 
because they make real time traffic simulations possible (see [|IJ and references therein). 
The model studied here is one such; for example, if the injection probability a is large the 
model may be regarded as describing the situation, familiar from everyday experience, of 
the reduction of a two-lane to a one-lane road by, e.g., the presence of construction work 
on one lane. It is a special case of the well-known Nagel-Schreckenberg model, obtained 
by requiring that the parameter t> max of that model satisfy t> max = 1 so that cars move at 
most one lattice spacing at each integer time. In realistic computer simulations of highway 
traffic, the Nagel-Schreckenberg model is usually used with v m3X = 5. However, in the case 
of OBC the phase diagram and density profiles are essentially independent of t> max |§, and 
in general it has been observed that, for modeling city traffic, it is sufficient to set t> max = 1 

The model is in fact a synchronous update version of the asymmetric exclusion process 
(ASEP), widely studied in both the physics and mathematics literature || || [7|] . The ASEP 
was originally introduced as an interacting particle system evolving in continuous time; this 
evolution is equivalent to the random sequential update (RSU) procedure, in which randomly 
chosen particles hop one at a time. Other updating schemes have also been introduced, 



including sublattice-parallel [§, U and ordered sequential procedures |TtJ, and the fully 
parallel updating (PU) scheme, which corresponds to the probabilistic cellular automaton 
described above. See Section 9 and for precise definitions and a review of current 



knowledge about these models. The ASEP with RSU and OBC has been exactly solved 
fT2| , |13| , 14] , using (among other methods) the matrix product ansatz, and these results have 



been extended to the sublattice-parallel and ordered sequential updating schemes [15, 10, 16 



II]. The model with PU has proved to be less tractable; for example, parallel updating can 
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induce strong short range correlations [|17| (these are absent under other updating schemes). 
On the other hand, parallel updating is important in practice for traffic modeling, both for 
efficiency — the PU is usually much faster than the RSU — and for effectiveness; for example, 
the Nagel-Schreckenberg model is always implemented with PU, since that method has been 
found to give the best agreement with measurements on freeway traffic p|. 

We remark that if one writes (3 = p/3 and a = pa and then takes the limit p — > 0, the 
model reduces to the ASEP with random sequential updating and injection and extraction 

rates (3 and a, respectively. 

Let us now briefly discuss known results which are related to our model. For the parallel 
dynamics on a ring (that is, with periodic boundary conditions), the exact solution was 
given in |2| ; here, in contrast to other updating schemes, the stationary state is not a simple 
product measure — occupation numbers at distinct sites are correlated. For example, if p 
is 1 and the density is 1/2 then the stationary state consists of free flowing particle-hole 
pairs, i.e., there is a strong particle- hole attraction. For p = 1, the dynamics is equivalent 
to rule 184 for cellular automata, for which transient properties have been analyzed 



The steady state for arbitrary p and overall densities is obtained by factorizing the weights 
of the configurations into clusters of length two the strong short range correlations 



persist. The steady state for the generalization of the model where each particle has its own 



hopping probability has also been solved JT9] 



Tilstra and Ernst ]^U] studied the case of OBC and p — 1. They obtained results which 
they argue to be asymptotically (i.e., in the limit of large system size) correct. In [11], the 



system was found to be exactly solvable on a special line in the phase diagram; from this 
special case and extensive Monte Carlo simulations, the phase diagram and formulae for the 
current and the bulk densities were conjectured. 

We now discuss briefly the nature of our solution. For the random sequential model, 
the initial breakthrough was the observation that there exists a recursion relation relating 
steady state weights (unnormalized probabilities) for a system of size N to those for a system 



of size N — 1 [12]. Equivalently, one may write the weights as matrix elements or traces 



of products of operators; requiring these operators to satisfy certain algebraic rules then 



implies that the weights satisfy the recursion relations [13]. The matrix product allows a 
more direct calculation of steady state correlation functions than the recursion relations 
ra, £TL HI. 



For the present model we have followed a similar line of attack. When we write the 
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weights as operator products, however, we must require that the operators satisfy quartic 
algebraic relations, which relate a product of four operators to sums of products of three 
and two operators (see section 2). This is in contrast with the quadratic relations found in 



previous works |13], [22|, ^3], |24|, |5|, |26], |27|, |29], |30fl . For example, the recursion relation 
for the weight f^(. . . 0100 . . .) relates systems of size N to systems of size N — 1 and N — 2 
in the following way: 

f N {. . . 0100 . . .) = (1 - p)/v-i(. . . 010 . . .) + f N -i(- . . 000 . . .) + pSn-2{- . . 00 . . .) . (1.1) 

We believe that the above method (recursion relations in the system size) should be of 
general interest as an analytic approach to probabilistic cellular automata [|3l], 32, 33], for 
which there are notoriously few exact results (see Schadschneider in [|TJ). 

We now summarize the content of the paper. In section |2| we define the model and 
provide the algebraic rules for a matrix product solution. The next section gives the proof 
that they indeed describe the stationary state; the argument proceeds by looking at blocks 
of consecutive particles and holes. In section ^ we show that the quartic algebraic rules may 
be reduced to quadratic rules by assuming the operators are two by two matrices whose 
elements are matrices, generally of infinite dimension, i.e., that the operators are rank four 
tensors, 

The reduction to quadratic algebraic rules allows us to relate the parallel update model 
to the model with other discrete-time updating schemes. Specifically in section 5 we show 
that the current and density profile for parallel update are simply related to those quantities 
for ordered sequential and sublattice parallel updating, although the relation between higher 
order correlation functions is more complicated. Thus, in solving exactly the parallel model 



in sections M — 10, we also obtain new exact results and prove conjectures for the other 



discrete-time models, for which only the asymptotic current was previously known. 

In section |6| several explicit representations of the matrices are constructed. As a first 
application, we solve the case p = 1 by means of 4 x 4 matrices in section ^. A detailed 
analysis is made of the two point correlation functions in order to highlight the oscillating 
decay of the correlation function which is a particular feature under parallel dynamics. 

We then turn to the task of obtaining the exact solution for general p < 1. The current 
phase diagram is derived using generating function techniques section R, and the asymptotic 
behavior of the density profiles in section |9] again using generating function techniques. The 
relevant Tauberian theorems are presented in appendix A. For finite systems, we calculate 
exact combinatorial expressions for the density profiles and two point functions (section [II]) . 
Technical details of the computations are contained in appendix B. By combining the results 
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from the two preceding sections, we compute the asymptotic bulk densities in section |TT]. A 
discussion closes the paper. 



2 Model Definition and Steady State Recursion Rela- 
tions 

In this paper we study the asymmetric exclusion process with parallel dynamics and open 
boundary conditions. We consider a one dimensional lattice, with N sites labelled 1 through 
N. Each site i may be occupied by a particle, in which case a binary variable r, satisfies t% = 
1, or empty, in which case T{ = 0. The n-tuple r = (r 1; . . . , t n ) specifies the configuration 
of the system. The dynamics is defined by requiring that at each time step three things 
happen: (i) all particles on sites 1, . . . , N— 1 with an empty site in front of them attempt to 
hop forwards, succeeding with probability p; (ii) if site 1 is empty then a particle attempts 
to enter the lattice there, succeeding with probability a; and (iii) if site iV is occupied then 
the particle there attempts to exit the lattice, succeeding with probability (3. All of these 
processes are stochastically independent. Note the particle- hole symmetry: the removal of a 
particle at the right end can be viewed as an injection of a hole, so the dynamics is invariant 
under the combined operations of interchange of i and N — i + 1, interchange of particles 
and holes, and interchange of a and f3. For example, we have 

(rj) N (a,(3,p) = 1 - (T N+1 _i) N (p,a,p) . (2.1) 



For p, a, and (3 nonzero, the configuration 1010 ... can be reached from any other, 
so the model, viewed as a finite state Markov chain, has a single irreducible component 



and hence a unique steady state |3J, which we denote by Pn( t ) is the probability of 
finding a system of size N in configuration r in the long time limit. In calculating Pn{t) it 



is convenient, as noted in earlier work on the random sequential model |12], [H], [13], |14|, first 
to define unnormalized weights /jv(t) and then to recover the probabilities via 



Pn{t) = f N {r)/Z N 



where 



z n = X1/jv(t) 

T 

the sum taken over all configurations of size N. 



(2.2) 



(2.3) 
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The idea is now to introduce a matrix product ansatz by writing 



f N (r) = (W\([[(l - Tj)E + Ti D)\V) . (2.4) 



This is to be read as a product of operators E and D (an E for each empty site and a D for 
each occupied site) contracted with a vector \ V) and dual vector (W\, yielding a scalar steady 
state weight; for example, ifiV = 6andr = 010001 then f N (r) = (W\EDEEED\V). 
This method originated in work on the random sequential model in which the operators 
were represented by infinite dimensional matrices. Later papers generalizing the original idea 



have employed operators represented by finite dimensional matrices [ 35] , j36|j or higher rank 



infinite dimensional tensors fl37| , |38f ; both these approaches will be used here. 

The operators E and D and vectors | V) and (W\ are required to satisfy certain algebraic 
rules, listed below. We determined these rules by finding the steady state explicitly for small 
system sizes and then guessing. We prove in the next section that these rules do indeed 
imply that the steady state of the system is given by fl2~2] ) and ( |2.4| ) for general N, and in 
section |6| we construct an explicit representation, thus verifying that the rules are consistent. 

The rules for the bulk are 

EDEE = (1 — p)EDE + EEE + pEE , (2.5) 

EDED = EDD + EED + pED , (2.6) 

DDEE = (1 - p)DDE + (1 - p)DEE + p(l - p)DE , (2.7) 

DEED = DDD + (1 — p)DED + pDD . (2.8) 

We also have rules involving three sites next to each boundary, 

DDE\V) = (l-/3)DD\V) + (l-p)DE\V)+p(l-P)D\V), (2.9) 

EDE\V) = (1 — (3)ED\V) + EE\V) + pE\V) , (2.10) 

(W\DEE = (l-a)(W\EE+(l-p)(W\DE + p(l-a)(W\E, (2.11) 

(W\DED = (1 -a)(W\ED+ (W\DD+p(W\D, (2.12) 

and two sites next to each boundary, 

DD\V) = P ^~^ D\V), (2.13) 
ED\V) = ^E\V), (2.14) 



(W\EE = ^£ — ^-(W\E, (2.15) 
a 

(W\ED = -(W\D. (2.16) 

a 

These rules permit the computation of all /at(t) (up to an overall constant, which must be 
assumed to be nonzero). However, a rather indirect argument, which we now discuss, is 
required when iV = 1 or iV = 2. It is convenient, and represents no loss of generality, to 
assume that 

(W\V}>0. (2.17) 

We may simplify (W\ED\V) using either ( p.!4j ) or ( 2.14Q ; equating the results shows that 
a{W\E\V) = (3{W\D\V), so that we may write 

(W\D\V) = £y{W\V), (2.18) 

(W\E\V) = -i(W\V), (2.19) 
a 



for some constant 7. Similarly, (W\DED\V) can be simplified using either ( [2.12 ) or ( |2.14| ), 
and this leads to 

(W\DE\V) = (1-P)(W\D\V) + (1 -a)(W\E\V)+py(W\V) . (2.20) 

Relations (|2.5|) - (|2.20|) allow the straightforward computation of all /at(t)- 

Equation ( |2.5| ), when inserted in ( j2.4| ), leads to the recursion relation ( |1 . 1| ) . From the 
set of all the algebraic rules one may similarly construct a whole set of recursion relations 
which uniquely specifies the steady state weights. It is more convenient, however, to work 
directly with the operator product. 

The algebra is not well defined when a or f3 vanishes, and the model is not interesting 
when p = 0. However, one may consider the random sequential limit discussed in the 

introduction, a = pa, j3 = p/3, p — > 0. Assuming that the operators E, D and vectors {W\, 
\V) of our algebra have limits E, D, (W\ and \V) under this scaling, and also assuming 
that 7 has the limit 7 = 1 (as it is true for the representations we construct in section |5|) 
one finds that ( |2.5| - |2.l6D and ( |2.18| - |2~2"0"D are, in the limit, consequences of the quadratic 



algebra of [0 , which reads 

DE = D + E, (2.21) 



(W\E = -(W\ 
a 



(2.23) 



Remark 2.1: The relations ( |2.5| )-( |2T8| ) can be used to obtain the steady state for our 
model with periodic boundary conditions (i.e., on a ring) if ( ^.4[ ) is replaced by /at(t) = 

Trf(n£i(l — Ti)E + TiD)^j . However, the algebra is not needed in this simple case, because 

the steady state is already known p. 



3 Proof of stationarity 



In this section we show that the operator algebra ( |2.5| )-( |2~T2] ) may be used to compute the 
stationary state of the ASEP with parallel dynamics. An elementary recursive argument 
shows that the weights /jv( r ) defined by ( [2.4p and the constant 7 introduced in ( [2.18| ) and 



fl2.19| ) satisfy /at(t)/7 > for iV > 1, so that ( |2.2j ) defines a probability distribution on the 



set of all system configurations. We must show that this distribution is invariant under the 
dynamics. To avoid consideration of many special cases it is convenient to first rewrite the 
algebraic relations satisfied by D, E, \V), and (W\ in more unified form. 



Note first that relations Q-(]2]i) , (^9|)-( ^l2D , and (|2~20l) all involve four factors (from 



among D, E, (W\, and \V)) on their left hand sides. These relations may be expressed by 
the single equation 

XDEY = a(XDY) XDY + a(XEY) XEY + p a(XY) XY . (3.1) 

Here X denotes either {W\, D, or E and Y denotes D, E, or \V). The coefficient a(S) of a 
term S (S = XDY, XEY, or XY) on the right hand side of ( |3.1| ) is determined only by S 
itself, not the relation under consideration: 



' l- 


P, 


if S contains DE, 


1 - 


a, 


if S contains (W\E, 


1 - 


/3, 


if S contains D\V), 


7, 




if S is (W\V), 


I !. 




otherwise. 



(Equation ( |3.2j ) represents a slight abuse of notation, since a(S) really depends on the 
form of S rather than on the value of S, which may be an operator, a vector, or a scalar; 
no confusion should arise.) We obtain another form of the relation (|3.1|) by first writing 
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XDEY = (1 - p)XDEY + pXDEY and and then using (|OD in the second term: 
XDEY = a(XDEY)XDEY+pa(XDY)XDY+pa(XEY)XEY+p 2 a(XY)XY . (3.3) 



The relations Q2.13 )— ( 2.16Q and ( 2.18|) — ( 2.19 ) can be similarly unified: 



XD\V) = £a(X\V))X\V), (3.4) 

(W\EY = -a((W\Y) (W\Y. (3.5) 
a 

Here, as in (|3.1| ), X is (W\, D, or E and Y is D, E, or | V). A second form of the relations 
and is obtained as was 0, starting from XFF = (1 - [3)XFY + /?XFF for 
and XFY = (1 - a)XFF + aXFF for 0: 



XD|^) = a(XD\V))XD\V)+pa(X\V))X\V) , (3.6) 
(W|£?y = a((W\EY) (W\EY + pa((W\Y) (W\Y . (3.7) 



Remark 3.1: (i) The set of terms on the right side of (|3.1|) is obtained by omitting either or 
both of the middle two factors in XDEY , and the set on the right side of (|3.3j ) by omitting 



neither, either, or both. Similarly, terms on the right side of ( |3.4|) and (3J3) arise through 



the omission of one operator, and those in fl3.6|) and (|3.7|) through the omission of zero or 



one. (ii) In ( |3.3| ), (|3.6|) , and (|3.7|) the power of p in each term is the number of operators 



omitted in obtaining that term; in ( |3.1| ), it is one less than that number. 

We now write down the stationarity condition which the weights /zv(t) must satisfy. If 
we imagine for the moment that our lattice contains two extra boundary sites, on the left 
and N + 1 on the right, then there are N + l bonds (i,i + l) across which an exchange might 
occur during one step in the evolution; here by an "exchange" across (0, 1) or (N, N + l) 
we mean the entry or exit, respectively, of a particle. Given a fixed configuration r, let us 
write A(t) for the subset of these bonds across which an exchange can occur in r and B{t) 
for the complementary subset of bonds across which an exchange might have occurred in 
arriving at r from some immediate predecessor; the bonds in A(t) correspond to (W\E, 
DE, or D\V) in the formula ( |2.4| ) for /jv(t), while those in £>(r) correspond to (W|.D, 
ED, or E\V). For C C B{r) write r c for the configuration obtained from r by making the 
exchanges corresponding to the bonds in C; the configurations r c , C C B(t), comprise all 
possible immediate predecessors of r. Then the stationarity condition has the form 

f N (r)= £ <C)f N {r c ). (3.8) 

ccb(t) 
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Here for any subset C of B(r), 7r(C) is the probability that, in the configuration r c , precisely 
the set of exchanges C (out of the set A{t c ) of possible exchanges) were in fact made; thus 
7r(C) is a product of the following factors: 

• a, if (0, 1) e C; (1 - a), if (0, 1) G A(t c ) \ C; 

• p, for each (i, i + 1) e C; (1 - p), for each (i, i + 1) e -4(r c ) \ C (1 < i < iV — 1); 

• /?, if (iV, iV + 1) e C- (1 - if (N, N + l)e A(r c ) \ C; 

with ^4(r c ) \ C denoting the set of bonds belonging to A{t c ) but not to C. 

To verify (|3.8|) we will use the relations above to transform each side to a common form. 
We need one more concept. The configuration r may be divided into blocks of successive 
zeros and ones; let S{r) denote the set of such blocks and for T C £(r) let tj? denote the 
configuration obtained from r by omitting one operator from each block in T . 

We first consider the left hand side of fl3.8|) . In the expression (|2.4|) for /at(t) we apply 
(p|) to each factor DE, to each D\V), and Q to each (W\E, if these occur. The 
resulting sum over all ways of omitting zero, one, or two operators at each of these bonds 
(see Remark 3.1.i) is equivalent to a sum over all ways of omitting zero or one operators 
from each block in r, so that from (|3.2|) and Remark 3.1 .ii we obtain 

£ 2^1(1 - af^{\ - p)v&(l - f3y^)f N _ m M, (3.9) 

where \T\ is the number of elements in T and x{T\ y{T\ and z(T) count respectively the 
number of factors (W\E, DE, and D\V) in the expression (|2.4j) for f{rjr). In the special 
case fN-\r\{ T F) — (W|l/"), which can occur only if r = 101010... or r = 010101..., 
there will be an additional factor of 7. It is important here that the coefficients a(S) in 
(P-3|), ( |3.6|) , and (|3.7| ) depend only on S, so that the coefficients in ( |3.9|) depend only on T 
and in particular are independent of the order in which the relations (|3.3|), ( |3.6| ), and ( |3.7| ) 
are applied at different bonds. Similar comments apply to the expansions using ( |3.1| ), (|3T^), 
and ( |3.7| ) which we will perform below. 

Consider now a weight /v(t c ) which occurs in the sum on the right hand side of (|3.8|) . 
Each bond in C corresponds in the expression ( p.4|) for this weight to a factor DE, D\V), 
or (W\E; we apply ( |3.1| ), ( |3.4| ), or ( |3.5|) to these factors. The result will be of the form 

f N {r c ) = Y,P m ^KF)fN-\r\{rr). (3.10) 
Since each of the relations we are using involves the omission of one or two operators (see 
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Remark 3.1. i), the JF occurring in ( |3.10|) will be those which, for each bond in C, contain 



one or both of the blocks abutting on this bond. The total number of operators omitted is 
\F\ and the number of times one of the relations was used is \C\, so that the factor pl^H^I 
is obtained directly from Remark 3.1.ii. The coefficient A(JF) is a product which contains a 
factor p/f3 if flS!]) was used (i.e., if (0, 1) G C) and p/a if (|3T5|) was used ((N, N+l))eC). It 
also contains factors arising from ( |3.2|) : a factor of 1 —p, 1 —j3, or 1—a for each application of 
fl3.ip, (p.4|), or ( p.5| ) which preserves or generates a factor DE, D\V), or {W\E respectively, 



and a factor 7 if fN~\r\{. T r) — (W\V). 

When ( ft.lOQ is inserted into the right side of (|3.8|) the double sum over C and T becomes 
a sum over all subsets T C S(t), 

Y: 7r(C)p^-M\(F)f N ^M, (3.11) 

F<ZS{t) 

since, given any J 7 , one may identify the corresponding C in (|3.10|) as the set of bonds in 
B(t) for which no block belonging to T abuts on C. It is straightforward to complete the 
proof by verifying that p^^^n^X^) is precisely the coefficient of fN-\r\{ T T) m (|3.9|) , 
and hence that ( |3.9| ) and ( [3.1 1|) agree. In particular, 7r(C)A(JF) contains a factor pl c l, that 
is, one factor of p for each bond in C; for internal bonds these factors are present in 7r(C), 
while if (0, 1) G C then vr(C) contains a factor a and A(JF) a factor p/a (the argument for 
(N, N + l) G C is similar). If (|3.9| ) contains a factor (1 — a), that is, if f{rj?) contains (W\E, 
then 7r(C) contains this factor if (0, 1) C and A(JF) if (0, 1) G C; the factors of (1 — p) and 
(1 — (3) in ( p.9|) are accounted for similarly. 



4 Reduction to a Quadratic Algebra 



In this section we show that the quartic algebraic rules (|2.5| )- (|2.16| ) can be reduced to 



quadratic rules by making a convenient choice for the operators involved. The trick is to 
write 

"-(£?)■ *-(«?)• ^ 

where £)i, D2, E±, and ^2 are matrices of arbitrary (in general infinite) dimension; that is, 
D and E are written as rank four tensors with two indices of (possibly) infinite dimension 
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and the other two indices of dimension two. Correspondingly, we write (W\ and \V) in the 
form 



(W\ = ((W 1 \, (W 2 \), \W) 



\Vx) 

\v*) 



(4.2) 



where (Wi\, (W 2 \, |Vi), and | V2) are vectors of the same dimension as Di and E\. We will 
show that the operators and vectors so defined satisfy the algebra of section ^| if D±, E\, 
(Wi\, and |Vi) satisfy the quadratic relations 



D 1 E 1 = (l-p)[D 1 + E 1 +p] , 



and D 2 , E 2 , (W 2 \, and \V 2 ) satisfy 

E 2 D 2 =p[D 1 + E 1 +p], 
E 2 \V 2 )=p\V x ) , (W 2 \D 2 = (Wx\p. 



a 



(4.3) 
(4.4) 

(4.5) 
(4.6) 



We now verify that imply (|2T5|)- (|2~20|) . First, by substituting Q) into the 

bulk relations (|2.5|) - (|2~lj|) one finds that to satisfy the latter equations it is sufficient that 



ExDxEt + E 2 D 2 E X = (1 - p)E 1 D 1 + (1 - p)E 2 D 2 + E\E\ + pE l , (4.7) 
E X D X E X D X + E 2 D 2 E X D X + E X D X E 2 D 2 + E 2 D 2 E 2 D 2 (4.8) 

= E 1 D 1 D 1 + E 2 D 2 D X + E 1 E 1 D 1 + E X E 2 D 2 + pEiDi + pE 2 D 2 , 
= {l-p)E 1 + {l-p)D 1 +p{l-p) , (4.9) 
D\E\D\ + D\E 2 D 2 = (l-p)E 1 D 1 + (l-p)E 2 D 2 + D 1 D 1 +pD 1 . (4.10) 

These relations follow from ([4.3|) and ( |4.5| ). For example, the left hand side of ( |4.7| ) becomes 

(l-p)E 1 [D l +E l +p)+p[D 1 +E l +p)E 1 = (l-p)E l D 1 +pD l E 1 + E 1 E 1 +pE 1 , (4.11) 



which another use of (|4.3j ) and ( [4.5|) shows to be equal to the right hand side. Similarly, 
with ( |4. 1|) and ( |4.2j ), the relations ( |2. 13| ) — ( |2 . 1 6| ) involving two sites next to the right hand 
boundary follow from 



D 1 D 1 \V 1 ) 
D 2 D 1 \V 1 ) 



P(H) ( 



£i£>i|Vi) + £7 2 £> 2 |yi) = ^Ei\Vx) + ^E 2 \V 2 ) 



(4.12) 
(4.13) 
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and these equations are in turn implied by Q) and O- The conditions (U) and (|27TO| ), 



involving three sites next to the right boundary, are obtained similarly, although the manip- 
ulations involved become rather tedious. Relations at the left boundary are obtained from 
symmetry considerations, completing the verification. 

In the remainder of the paper, we will consider only representations of the quadratic 
algebra having the form of ( f4.1|) and ( |4.2|) . As we now discuss, computations in this repre- 
sentation are simplified by the fact that all quantities of physical interest may be expressed 
in terms of D±, E\, (W\\, and \Vi). This also leads to connections with other updating 
procedures for the ASEP (see section ||). 

Let us first derive such an expression for the normalization constant Z^, given by 

Z N = (W\C N \V) , (4.14) 

where 

C = D + E=(% E Q 2 ) (4.15) 



and Ci = D 1 + E x . Now, 



r N_( G{N) G(N-1)E 2 \ 

-\D 2 G(N-1) D 2 G{N-2)E 2 )> l4 ' iDJ 

where by convention G(—l) = and 

N 

G(N) = ]T K N ~ n {-p) n (4.17) 



ra=0 



for iV > 1, with 

K = (d+p). (4.18) 

This may be proven by first checking the case N = (G(0) = 1) and then verifying the 
recursion G(N + 1) = dG(N) + E 2 D 2 G(N -1) = (K — p)G(N) +pKG(N — 1). Note that 

G(N) + P G{N - 1) = K N . (4.19) 

From (|4.16|) , ( |4.19| ), and the action ( |4.6| ) of E 2 and D 2 on the boundary vectors, we have 



(W\C n = ( (W^K™, (W 1 \K n ~ l E 2 ) , (4.20) 
K n \Vi) 



C n \V) = ( n .^K 1 l > T/ v ) , (4-21) 



13 



which leads to 



where 



Z N = Z N + pZN-l , 

z n = (W 1 \K n \V L ). 



(4.22) 
(4.23) 



An important quantity in determining the phase diagram is the current Jn, which is the 
probability that a particle passes through a particular bond in a particular time step. It is 
given by any one of three equivalent expressions: 



J N = - n)} = P<Ti(l - r i+1 )) = P(t n ), 
where 1 < % < N; the second of these may be written using the algebra as 



(4.24) 



Jn =P 



(WIC^DEC"- 1 - 1 ^) 



(4.25) 



Now (|4~20D and (pL2lD yield 



(W\C n D = ( (W 1 \K nr - 1 [D 1 +p], 
EC n \V) = 



[E 1 +p]K n - 1 \V 1 ) 




(4.26) 
(4.27) 



and the algebraic rule (|4.3| ) implies that 



[D 1 +p}[E 1 +p]=K, 



(4.28) 



so that from ( [4^5] ), ( ggg ), and ( gg3p 



J. 



PZN-1 



iV 



(4.29) 



This expression again involves only the matrices E\ and D\ and vectors and |Vi). 
We may similarly express the one-point correlation function or density profile, 



(r l ) N = ^(W\C l - 1 DC N - l \V), 



(4.30) 



in terms of E\, Di, (Wi\, and |Vi): 



(r<> 



A' 



{W^K^jDi+pjK^lVj) 
z N +pz N _i 
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(4.31) 



where we have employed ( |4.26|) and ([4.21 ). Similar expressions are possible for higher 
correlation functions; for example, the two-point correlation function (rj(l — Tj))n is given, 
from (p^D , and ( ^27|) by 



JV 



-^(W 1 \K i ' 1 (D 1 +p)G(j - ^ - +p)K Jv ^'|y 1 ) . (4.32) 



Because G(n) is an alternating sum, ( f4.32| ) is more complicated than the corresponding 
expression in the random sequential case; this reflects the fact that stronger correlations 
exist under parallel dynamics than under random sequential dynamics. 



In section 10 we will give a more explicit formulae for taW; see (|10.12|). 



5 Relations with other models 

The reduction to a quadratic algebra and the expressions for the current and correlation 
functions, derived in section |], lead to relations between the ASEP with parallel updating 
and the same model with certain other discrete-time updating procedures. The procedures 
in question can in fact be defined for more general site variables and for any local dynamical 
rules which assign to each configuration of a pair of sites at time t a new configuration at 
time t + 1 with some given probability, and similarly to each configuration on the leftmost 



or rightmost site a new configuration. We recall these procedures briefly; see [[Tl|] for precise 
definitions. In the ordered sequential update sites are updated one at a time, starting the 
right end of the system and proceeding sequentially to the left end (backward ordered), or 
vice versa (forward ordered). In the sublattice parallel update, all site pairs i,i + 1 with i 
even are updated at one time step and and all such pairs with i odd at the next time step. 
Let us denote these procedures by the symbols T<_, TU and T sp (more precisely T<_, Tl, and 
T sp are the transfer matrices for the different procedures, however since in this discussion we 
do not require any properties of the transfer matrices, we do not give detailed definitions). 



It can be shown ]39] using the matrix product formalism that in general the procedures 
T<_, TU and T sp lead to stationary states which may be regarded as physically equivalent. In 

particular, the current is independent of the update procedure; we write J # for this common 
value: 

Jn = Jn = Jn = Jn ■ (5-1) 
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The density profiles are also closely related (we use the notation of the ASEP, but a corre- 
sponding result holds in general |39fl ): 



/ \ S p _ J ( T i)N, i even, /^\sp*_J ( r «)w") * odd, , > 

{Ti)n ~ \ (rfc, i odd, {Ti)N ~ \ ( n )j;, i even, 

where (jif^* is the density for T sp after the first (even) sublattice has been updated. Similar 
relations hold for some higher order correlation functions. 

More is known in the special case of the ASEP, to which we limit our discussion in what 



follows, based on simpler realizations of the matrix product ansatz in that case [151, [10], [16 



TTj| . In particular there is another relation between density profiles, 

( T i)N = ( T i)N ~ (5-3) 

which, with ( |5.2|) , shows that any one of (tj)^, (Ti)]y, and (rj)^ determines the others. 
Moreover, the asymptotic current liniAr^oo J* (and therefore the phase diagram) is known 
fl6fl, and so are the density profiles and correlation functions in the case p = 1 |L5], 10]. 



A fully parallel updating procedure is not naturally defined for arbitrary local dynamical 
rules, due to the possibility of conflict when the rules are applied simultaneously to pairs 
overlapping of sites, but for the ASEP such a fully parallel procedure, which we will denote 
by T\\ , has been defined in section |2[. We now want to relate this model to the ASEP with 
update procedures T_>, T<_, and T sp , considered above; for the ASEP with parallel update 
we work with the reduced version of the algebra established in section |], in which everything 
is expressed in terms of the matrices E\ and D\ and vectors (Wi| and |Vi), which satisfy 
the algebraic relations ( |4.3| , f4.4| ). The idea is to introduce new operators e and d by 



d = D 1 , e = E 1 +p, (5.4) 



so that K = Ei + D\ + p = e + d. From (|4.3|, |4.4j), the new matrices satisfy 



de = d+(l-p)e (5.5) 

and 

d\V 1 ) = p ( 1 ~® \V 1 ) , (W 1 \e=(W l \^- . (5.6) 
p a 



Surprisingly, these equations are precisely the algebraic relations [|TU| for the matrix prod- 



uct solution of the ASEP with forward updating, that is, the steady state weight for the 
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configuration r with updating T_> is given by 

N 

(W 1 \([[(l-T^e + nd)\V 1 ) . (5.7) 
i=i 

Writing c = e + d we see that the normalizing factor for T_> is 

(W^IVJ = {W 1 \K N \V 1 ) = z N , (5.8) 
a constant already introduced in (|4.23|) , and the current is 

Z N Z N 

The density profile is given by 

. (H'- 1 |c i - 1 dc A '-'H' 1 > (W 1 \K i - 1 D l K"- i \V l ) ,, 

<n> = — t n — = — www, — • (5 0) 

These formulae imply certain simple relations between physical quantities in the parallel 
and the ordered sequential ASEP. From ( |5.9| ) and the formula (|4.29 ) for the current in the 
ASEP with update Ty, Jn = V z n-i/{ z n + P z n-i), we have 

Jn = -4L • (5.11) 

1 + J N 



In particular, Jn < «/*, a result which is intuitively clear. Similarly, comparing ( |5. 10 ) with 
the formula ( [4.31 ) for the density profile for the parallel ASEP leads to 



l + J* " l + J N 



T i)N = „ , T jt (5.12) 



where we have also used ( |5.3|) . It is also possible to derive similar formulae for the higher 
correlation functions although these are not so simple. For example, the two-point function 



( 1.32 ) for Tii can be written as an alternating sum over two-point functions of T_ 



1 

(T;(l - Tj)) N = — Yj (-P) n ( r i(! - T j-n))N-n Z N-n ~ Pi^N^^ZN-n-! , (5.13) 
Z N n=0 

although this formula appears to be too complicated to be of practical use in obtaining 



(r»(l — Tj)) N from (rj(l — Tj))^ or vice versa. We emphasize that the formulae ( |5.11| ), 
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( |5.12| ), and (|5.13|) hold for all a, /3, and p. Of course, using the formulae established earlier 
in this section, we obtain relations for the one- and two-point functions with updates T«_ 
and T sp . 

The relation ( |5.11| ) for the currents and the relation for the bulk densities which follows 
from ( |5.12| ) were already conjectured in |Tl[ (see Table 1 in that reference). 

These relations between different models have various consequences. On the one hand, 
established results for T«_, T_ and T sp serve as a check for the results we will derive later; 
this applies to the asymptotic values of the currents and hence to the phase diagram which 
we derive in section ||, and to some of the results in the special case p — 1. 

On the other hand, and more importantly, most of the results that we will derive in later 
sections — explicit representations of the algebra, some of the detailed results for the case 
p = 1, asymptotics of density profiles, and finite volume formulas for current and density 
profiles — apply to all of the update procedures Ty, T<_, T_», and T sp . We will state these 
in terms of Tji, since that procedure is the main focus of this paper, but results for other 
updates are easily obtained via the formulae of this section. 

Finally, we want to point out a further connection to another known model, the ASEP 
with random sequential update on a ring with a second class particle |[24, 25 1. In the ASEP 
with two species of particles a first class particle hops onto a vacant site to its right with 
rate one and exchanges positions with a second class particle to its right with rate r, and a 
second class particle hops onto a vacant site to its right with rate s. Let us summarize some 
known results for this model. The stationary state on a ring can be written as a matrix 
product state [p4| ; the corresponding matrix algebra is given by (|2.21|) -( [2.23|) , where now 

the matrix for first class particles is D, the matrix for holes is E, the matrix for the second 
class particle is |V r )(W / |, and 5 = r, (3 — s. Since we are considering a closed system, it is 
convenient to work in a grand canonical ensemble pi] with fugacity x for first class particles. 
When the system contains just one second class particle, the density of first class particles 
i sites ahead of the second class particle in a ring of iV sites is given by 

N (W\C N ~ l \V) ' v ' ; 

where C = xD + E, with x the fugacity. 

Now we can map the algebra ( |5.5|) , ( |5.6|) onto that of ( |2.2i|) -( p.23| ) by defining 

(l-p)D = d, E = e, (5.15) 
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and 



a 



P 

a 



(3 



13(1 -p) 



(5.16) 



p(l-(3)- 

Then c = e + d is equal to C\ x=1 , so that one can obtain the density profile for the ASEP 



with T_>, and hence for the ASEP with Ty, from the grand-canonical density profile ( p.!4| ) 
seen from the second class particle: comparing ( |5.10| ) and ( |5.14| ) yields 



(i-p)(n)N 



( T i)iV+l 



r=a/p, s=/3(l— p)/p(l— 0), x=l—p 



(5.17) 



Other quantities in the models can be related similarly. When r = s = 1, (r,)^ p is known 
exactly for both finite and infinite systems p4| , yielding (r^)^ for the case a = (3 = p. The 



single second-class particle model for general r and s has been studied |25j in the canonical 
ensemble of a fixed number of particles on the ring; to the extent to which the canonical 
and grand canonical ensembles are equivalent in the large N limit, these results should 
correspond with our results for the model with parallel update. We checked that indeed the 
phase diagram we will derive in section [| translates to the correct phase diagram for the 



model of 25 . 



6 Representations of the quadratic algebra 

In this section we discuss several explicit representations of the quadratic algebra (|4.3|) , 



(Q[). The situation is like that for other, similar matrix product algebras: finite dimensional 
representations exist for a few special parameter values, and some representation, typically 
infinite dimensional, exists for all values. 



In ]TI] it was observed that for parameter values on the line 

(l-a)(l-(3) = l-p, (6.1) 

the weight of a configuration in the stationary state could be written as a product of factors 
corresponding to clusters of length two. Thus, we expect a simplification of our algebraic 
relations (|Q|)-(fO|) along fl6.1|). Indeed, for this special case the matrices E\, D\, E 2 , D 2 
and the vectors (Wi\, (W2I, \ Vi), | V2) can be chosen to be scalars: 

^ 1 — ^ 1 — a . . 

Di=p— f- , E 1= p , (6.2) 

D 2 = E 2 =pv, (6.3) 
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\Vi) 

\v 2 ) 



(W 1 \ 
(W 2 \ 



1, 



(6.4) 
(6.5) 



where v = yp/aj3, so that E and D are 2x2 matrices. 



In section |7| we will treat the case p — 1 (with a and /3 arbitrary), in which one can 
choose a two dimensional representation of the Ex, Di algebra. Here, the bulk dynamics 



is completely deterministic. However, the physics is still interesting |[20|| , since a and (3 
can induce different phases. The algebra is sufficiently simple that we can derive explicit 
formulae for quantities such as the fluctuations in the number of particles in the finite 
system. 



Finite dimensional representations of the matrix product algebra ( |4.3| ), (fO|) exist only 
along the line ( |6.1| ) and when p = 1. This can be seen by using the mapping ( |5.15| ), ( |5.16| ) 
in section g, because it has been proven [|X3|] that the matrices in (|2.21| )- (|2.23|) have to be 

infinite-dimensional except in the case a + (3 = 1, which is by ( |5.16| ) equivalent to condition 
( 16.11 ). Note that the mapping ( p. 15) ) and ( |5.16| ) is not well-defined for p — 1; however, in 
this case we know that there is a finite dimensional representation, because we construct it 
(see section 0). 

Let us now turn to the case of general p,a,/3, in which the representations must be 
infinite dimensional. Such representations are of use both as a calculational tool and also 
to demonstrate that non-trivial representations of the algebra actually do exist. It can be 
shown by direct calculation that the following expressions satisfy 



(6.6) 



= ( 1, 0, . . ) , |K) = ( 1, 0, . . ) 



E 1 





b 










(1 


-p) {l-pfl 2 










(l_p) 










(1-p) . 




V • 






■ ) 


/ p{l — a) /a 








• \ 


b 


(1 


-p) 







(1 


_p)i/2 (!_ p ) o 










(l-p) 1/2 (l-p) 




V • 






• / 



(6.7) 



(6. 
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where 



p_ 

af3 



[(l-p)-(l-a)(l- 



(6.9) 



and 



E 2 = gD 1 , D 2 = gEt , 
\V 2 ) = /3/(g(l - P))\Vi) , (W 2 \ = (Wi|a/(^(1 - «)) 



(6.10) 
(6.11) 

where g = Jp/ (1 — p). These representations are not defined when p = 1, however, as 
explained above, we will treat this case in section [7]. Note that for (1 — a)(l — (3) = 1 — p 
we have b 2 = 0, so that the (1, 1) elements of the matrices decouple from the rest and we 
are left with 2x2 operators D,E (see also 4.2). This proves a conjecture of [ll|] for the 
structure of the steady state. In the limit p — > one recovers the representations (36) and 
(37) of M . 



7 The case p = 1 

In the case p = 1 the hopping of particles in the bulk is deterministic; the only source of 
randomness comes from the parameters a and (3. We shall see that a and (3 can induce 
a high density phase (a > (3) and a low density phase (/3 > a). For a = (3, these two 
phases coexist (see below). Note that the recursion relation following from the algebraic 
rule ( |2.7| ) implies that the weight for configurations containing a particle-particle-hole-hole 
string is exactly zero. This is also immediately evident from the dynamical rules, since these 
substrings can never be created (they are "Gardens of Eden" that can never be entered once 
left). 

One can choose 2x2 representations for the operators E\ and Di, so that E and D 
are 4x4 matrices. In particular, one can verify that the following explicit representations 
satisfy O-O: 

(Wx\ = ( 1, ) , \V X ) = ( 1, ) t , (7.1) 

Z^(« 7) , J) , ( 7 . 2 ) 

(W 2 \ = ( 1/a, ) , \V 2 ) = ( 1/a, f , (7.3) 
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where 



(l-g)(l-/3) 1 

c= y 3 ' a = \—s ■ ( 7 - 5 ) 



a/3 y a(3 



In order to compute expectation values, it is convenient to diagonalize K (see ( [4.18| )). 
It turns out that this can be done for a ^ j3; we discuss this case first, but omit details of 
the diagonalization. However, note that the eigenvalues of K are 1/a and 1/(3. 

Using the diagonalization, it is straightforward to compute = (Wi\K N \Vi) by writing 
|Vi) as a superposition of eigenvectors of K and one obtains 

Q.-P)al3- N -(l-a)l3a- s 
z N = — . (7.6) 



Equation (|4.22|) then leads to 



(1 - (3 2 )a (3~ N - (1 - a 2 )(3 



a: 



-N 



Thus, the current J N follows from ( |4.29| ) 



JN - ap (l-{p) a ™-(l- a *)0™ • (7 - 8) 

From this equation, it is clear that there are only two phases (for a ^ (3), a high density 
phase when a > (3 and a low density phase when (3 > a. As N —>■ oo, Jn approaches 
(3/(1 + (3) or a/(l + a), respectively. These results were conjectured in [[□], |2C| . 



With the diagonalization and ( |4.31| ) it is straightforward to work out a formula for the 
density profile (Ti) N (a,j3) (i = 1,...,N): 

a N +\l -0)- g*+iq(l - a) - (P/aY(l - a)(l - (3)a N+1 
{Tl)N (1 - /3 2 )a N + 1 - (1 - ot 2 )(3 N + 1 ' 1 ' 

It follows that in the limit iV — > oo the density profile in the high density phase a > (3 is 
given by 
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This implies that the density is constant for i ^> 1 (in particular, in the bulk and near 
the right boundary), and decays exponentially near the left boundary (on a scale given by 
log (a//?)) towards the bulk density: 

(r)^ = YTp ■ (7 ' n) 

The corresponding formulae for the the low density phase (3 > a can be obtained easily 
directly from ([7.9|) or via the particle hole symmetry Q2.1|) . One obtains 



(r> bu ik = — - , (7.12) 
1 + a 

and an exponentially growing density profile near the right boundary. These values for the 
bulk densities are in agreement with conjectures in |TT], 



The two-point correlation function (for a ^ (3) can be computed via ( [4.32 ) and is given 
by the following expression: 

( TiTj ) = f(a,^N)- 1 x { (3 N+1 a 2 {a-l){l + (3) + a N+1 {l + a){l-f3) (7.13) 

+P N+1 a(a - 1)(1 + P)(-a) j - 1 + a N+1 (l + a){j3 - 
+a N+1 {a - 1)(1 + a)(l - /3)(/3/a) i 
+a N+1 {a -1)(P- l)(a - (3){-(3y{-a)-' 1 
+a N+1 (a-l)(l-f3)(l + f3)a((3/a) j } , 

with i < j and 

/(«, P, N) = (1 + a)(l + p) [ (1 - /3 2 )a^ +1 - (1 - a 2 )^ +1 ] . (7.14) 

From this one can derive asymptotic expressions; for example, in the high density phase 
(a > (3) the bulk two-point correlation function (J = i + r, 1 << i < i + r) is 

{r t T i+r ) hulk = ( 1 + p\ 2 ■ ( 7A5 > 

Note the oscillating nature of ( |7.15|) , which can be interpreted as a particle- hole attrac- 
tion which is created by the simultaneous updating (see introduction). The corresponding 
truncated correlation function, 

g(i,j) = (nrj) - foXr,-) = , (7.16) 
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decays for (3 ^ 1 exponentially to zero on a scale 1/| ln/?|. 

Let us now turn to the case a = (3. In that case, K cannot be diagonalized, but there 
exists a similarity transformation which reduces K to Jordan normal form and makes all 
the necessary calculations straightforward. We just list the results: 



z N = ^[N(l-a) + l 



Zjsr — — - 



a 



N 



J iv = a 



N(l - a 2 ) + a 2 + 1 

N(l -a) + a 
N(l-a 2 ) + l + a 2 



and 



which yields for large iV 



(r*> 



N 



[1 - a) 2 i + aN(l - a) + a 
N(l -a 2 ) + 1 + a 2 



(7.17) 

(7.18) 
(7.19) 

(7.20) 



1 + a 



(7.21) 



The error term is small when N(l — a) ^> 1. Again, these expressions coincide with the 
expected formulae The linear profile in Q7.21]) can be interpreted as arising from a 

uniform superposition of states with localized shocks. 

The two point correlation function is now 

( TiTj ) = t(a, N)- 1 x {(1 - a)(l - a 2 )i + a(l - a)(l - a 2 )j + Na 2 (l - a 2 ) + 2a 2 

)(i - j) + Na(l - a 2 ) + a(l + a 2 ) ) } (7.22) 

where i, j = 1, 2, . . . , L, i < j and t(a, N) = (1 + a) 2 [N(l - a 2 ) + 1 + a 2 }. When a = 1 
this reduces to 

1 



( T i T j)a= 



l + (-l)^ 



(7.23) 



which is the expected result since then in the steady state only the two configurations in 
which there are no particle-particle or hole-hole pairs occur with non zero probability. The 
alternating structure of (|7.23| ) does not show up in the density profile (which is flat here) 
because these two configurations have equal weights. 
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For q^I, the truncated correlation function g(i,j) simplifies for large N: 



g{iJ) = jTTar^ /N){1 ~ U/N)){1 ~~ a)2 + ^ a)rail + (r/iV)) } + 0{l/N) ' (7 - 24) 

where r = j — i. Let us briefly discuss the behavior of g(i,j) for fixed i and j — i, . . . , N. 
The oscillating part of g(i,j), which is not present with other updating schemes, decays 
exponentially with j on a scale l/|lna|. Therefore, for sufficiently large j, g(i,j) decays 
linearly to zero with slope — (j^) 2 jf2- Figure 1 shows two examples of g(i,j) for a system 
of 100 sites. The strong oscillations present for a = (3 = 0.9 arise because the density in 
the system is nearly at its maximum value of 1/2, so that if a site is empty its two nearest 
neighbors are probably occupied. When a — (3 — 0.1, on the other hand, each site is for 
some typical configurations in a region of low density (if the shock is to its right) and for 
some in a region of high density (if the shock is to its left), so that the truncated correlations 
are positive. 

We now turn to the calculation of the fluctuations in the number M of particles in the 
system, still considering the case a = [3. We write 

N 

(M) = J» (7.25) 

i=l 

N 

A 2 = (M 2 )-(M> 2 = 2]T(T i T J > + (M)-(M) 2 (7.26) 

i<j 

and must sum up the expressions ( |7.20| ) and ( |7.22| ), respectively. The first summation is 
trivial: 

N 2 (l-a 2 ) + Ml + a 2 ) 



For a = (3 = 1 we obtain the expected result (M) a= i = N/2. The summation of ( [7.221 ) 



is more tedious. It is convenient to use J2fi<j h(r) = J2r=i(N — r)h(r), where h(r) is an 



arbitrary function of r. One is then left with well-known sums of the form J2^=i r k (— a) 
(k = 0, 1, 2). Altogether one obtains 
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A 2 (a,iV) = t(a,N)- l x{^ [1 + a - a* - a A ] - N 2 a 2 [- 



1 + a 



2Na 2 [ 2 + a + a2 ] + (l-a + a 2 )(—) 2 
1 + a 1 1+a 



— (—a 



2a " 2 {N(l-a 2 ) + l-a + a 2 ]} 



1 + a' 

+ (M) - (M) 2 . (7.28) 

For a = 1 ( [7.28Q yields 1/4 for iV odd, and for N even, as expected. 
We now take N ^> 1 and keep only the highest order in N. This gives 

12(l + aP ' 1 ' 

which can be rewritten, using (|7.21|) , as 

N 2 

A 2 = — (p right -p lcft ) 2 , (7.30) 

where p r ight(pieft) is the asymptotic density at position N(l) given by ( |7.21| ). This is precisely 
the result which is to be expected when one considers the linear profile as a superposition 
of uniformly distributed random shock positions (step functions). 



8 Derivation of the phase diagram for general p 

In this section we determine the asymptotic behavior, for all values of a, /3,p, of the quantity 
z N = EwvK N introduced in (|4.23|) and hence, through ( f4.29| ), of the current J^; the different 
possible asymptotic forms determine the distinct phases of the model. Our method is to 
study the generating function 

oo 

e (A) = £ X N z N . (8.1) 

N=0 



We will use the explicit representation (|6.6j )- (|6.8| ) of the operators D\ and E% and the 



vectors (Wi\ and |Vi), and will write |n), n = 0, 1, . . ., for the basis of the space on which 
D\ and E\ act, and (n\ for the dual basis, so that \Vi) = |0) and (Wi\ = (0|. Note that 
z = (W 1 \V 1 ) = l. 



26 



From ( |6.7|) and ( |6.8|) it follows that K = E\ + D\ + p has the form 



K 



( c b 

b 2-p (l-p) 1 ' 2 

(l-p) 1 / 2 2-p (l-p) 1 / 2 

(l-p) 1 ' 2 2-p 



\ ■ 



where c = p(a + (3 — af3)/a(3 and b is given by ( |6.9|) . From ( |8.2|) we find 



(<d\K N+1 \0) 
(l\K N+1 \$) 



and for n > 1, 



(n\K N+1 \o) 



c(0\K N \0) + b(l\K N \0), 

6(01^10) + (2-^(11^10) + (1 -p) 1/2 (2|i^|0> 



(l-p) 1/2 (n-l\K N \0) 

+ (2 - p)(n\K N \0) + (1 - pf/ 2 (n + 1\K N \0) 



If we now define the generating functions 



\ N (n\K N \0), 

N=0 

we easily obtain from (|Q|)-(|S75|), using (n\0) = 5 n o, that 



(l-cA)0 o = bXQx + 1, 
(l-(2-p)X)Q 1 = 6A6 o + (l-p) 1/2 A0 2 



and for n > 1, 



(1 - (2 - p)X)Q n = (1 - p) 1/2 Ae n _! + (1 - p) 1/2 A9„ +1 



The solution of (18.91) is 



6,, = Au n for n > 



where A is a constant to be determined from 08.81) and 



u 



l-\(2-p)- ^J(l + \p) 2 -4A 
2A(1 -p) 1 / 2 
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(the positive root is discarded being non-analytic at A = 0). Writing X = Au and 2 = Au 2 
in (|8.7|) and (|8.8|) , and eliminating A from the resulting equations, yields 



6r 



1 



1 -c\-b 2 Xu/{l -p)V 2 



(8.12) 



From ( |8.12| ), ( |8.11|) , and some tedious algebra, we find that 



e (A) 



a(3 [2(1 - p)(a(3 - p 2 X) - a(3b 2 (l - pX) - a(3b 2 ^J(l + pX) 2 - AX 



2p 4 (l-/?)(l-«)(A-A hd )(A-A ld ) 



(8.13) 



where 



Aid 

Ahd 



a(p — a) 
p 2 (l — a) 

P(p - g) 

P 2 {l-(3) 



(8.14) 
(8.15) 



Equation ( |3.13|) shows that Oq has square root singularities at the two points 



xL = (8 . 16) 

which, if we assume that the parameters a, /3, and p lie in the relevant range < a, /3,p < 1, 
are on the positive real axis. Thus Go is naturally double valued and has single valued 
determinations (branches) on each of two sheets — copies of the complex plane cut along the 
real axis between the two roots. We are primarily interested in the behavior on the first 



sheet, the plane on which, for A small and real, the square root in ( |3.13|) is positive and 
hence 0q(O) = 1 (see ( |3.1|) ). 6 also has two simple poles, at Ay and Am- 



As discussed in Appendix [A], the coefficients in the power series (|8.1| ) will grow as 
N 7 Xq , where A is the singularity of 0o(A) on the first sheet nearest to the origin (A > 1 
always) and 7 is determined by the nature of that singularity. Thus 



and hence 



lim = A , (8.17) 



lim J N = (8.18) 

AT-oc 1 + pX V ' 
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Three regions in the parameter space must be considered according to which of the three 
singularities is closest to the origin. As we shall see, the singularities A^ and Ahd may or 
may not be present on the first sheet of the complex plane. For the parameter values where 
one or both do occur they are closer to the origin than A~ c . It is convenient to introduce 
the quantity 



q = q(p) = -p (8.19) 

in discussing the resulting phase diagram (thus 2q — p = q 2 and A~ c = (q/p) 2 )- The phase 
diagram is shown in figure 2. 

(i) The maximum current region: q < a and q < (3 (region C in figure 2). 



For these parameter values the numerator in ( |3.13 ) vanishes at Ay and Ahd when the 



square root is positive: the poles lie on the second sheet and the singularity closest to the 
origin in the first sheet is A~ c . Then ( |A.7| ) implies that 

ZN = g V ™ i 0(N / (\~) ), (8.20) 



and hence from ( 8.1 



hm J N = PX ™ = WO . (8 . 21) 

*-oc N l+p\ mc 2 1 ' 



Note that the prefactor in fl8.20|) is singular at the boundaries a = q, (3 = q of the maximum 



current region; near these boundaries one needs larger values of N for the leading term in 
( |S.20| ) to well approximate z n . 

(ii) The low density region: a < (3 and a < q (region A in figure 2). 

In this region the pole Ay lies on the first sheet and is in fact the singularity of Go (A) 
closest to the origin, and from (|A.5|) , 

_ (3{p + a 2 -2a) 1 _ N 
[p -a){J3- a) (A w ) N 

for some s > Ay- Thus from ( |8.18| ), 

an j K = : ^L. = ^PZ^. ( 8 .23) 
n^oo 1 + p\ ld p — a 1 

(iii) The high density region: (3 < a and (3 < q (region B in figure 2). 
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The asymptotic behavior of zn here is obtained from ( |B.22| ) by interchanging a and (3 
and replacing Aid by Ahd- In particular, 



lim J N = ^> . (8.24) 

N-^oo p — fj z 



9 Asymptotics of the density for general p 

In this section we calculate the behavior of the particle density near the left end of the 
system in the limit of infinite system size and for all values of a, (3 and p; behavior near the 
right end can be recovered from the symmetry For m,n > let 

EwvK N ~ m ~ n D^ 1 K n , . 

t mjn = lim * 9.1 

and introduce the generating function 

*(x,j/)= E s'YW (9.2) 

m,n>0 



Our goal is to calculate ty x (0,y). For it follows from (|4.31| ) and (|8.17| ) that p n , the density 



at the (n + l) st site to the left of the right boundary in the infinite volume limit (where 
n — 0, 1, . . .), is given by 

y / \ v EwvK N - l - n (D 1+ p)K n ty n + P X 
p n = hm {T N _ n } N = hm ■ = — — — , (9.3) 

JV->oo N^oo Zn+PZn-1 1 + P^o 

so that the generating function $(?/) for the p n is 

Hv) = E y n Pn = — ^- E y n (ti,n + P \o) = — ^- f *.(o, v ) + . (9.4) 



Now t ,n = 1 for all ra, so that V{0,y) = (1 - and from (g7g ) and (|8.17|) , t m>0 



(X p(l-(3)/(3) m for all m, so that *(x,0) = P/(P-xXop(l-0)). From §J) it follows that 
_Di_K" = + (1 —p)K +pDi and this, together with (|8.17|) , implies that for m, n > 1, i OT)n 
satisfies the recursion 

tm,n ^m+1, n-1 + A Q (1 - p)£ m _i,n + A pt m ,„-i. (9.5) 
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Thus 

v(x,y) = ^ x m t mi0 + ^ y%,„ + E x m y n tm, n 

m>0 n>l m,n>l 

= 0) + [tf (0, y) - 1] + 2 + A (l - p)tm-l,r. + A pt m ,„-i] 



m,n>l 



= 0) + [tf (0, y) - 1] + y) - (0, y) - x* x (0, y)\ 

+A (1 - p)a:[tf y) - *{x, 0)] + X oP y y) - (0, y)]. (9.6) 

Multiplying (|9.6|) by — x, collecting terms, and using the relation \I/(0, y) — 1 = y\P(0, y), we 
see that \l/ satisfies the equation 

D{x, y)*{x, y) = A(x, y) + xy$ x (0, y), (9.7) 

where 

D(x,y) = \ {l - p)x 2 -{1- X oP y)x + y, (9.8) 

A(x,y) = -y[x(l-\ o p)-l]y(0,y)-x[l-\ o (l-p)x]y(x,0). (9.9) 

Now the branch of the curve D(x, y) = given by x = £_(y), where 



l-A oP y±^A(y) 

^ = 2A (l-p) (9 - 10) 

with 

A(y) = (l + \ py) 2 -4\ y, (9.11) 

are the roots in x of D(x,y) = 0, passes through the origin. But ty(x,y) is analytic at the 
origin, so that ( ^.7[ ) can hold only if the right hand side vanishes on this curve. This yields 
the desired equation for ty x (0,y): 

^-(y),y) _ e-(y)(l-Aop)-l /?[l-A (l-p)e_(y)] rqi9 s 

x[,y) Uv)v C-(y)(i-y) y[/3-^-(y)A p(i -/?)]" 1 j 



If we insert (|9.12|) into from (9^4) and rationalize the resulting expression we obtain 

$(y) 



2y - 1 + A py V A(y) 



1+pAol 2y(l-y) 2y(l - y) 



P{j>-P)yf*ti) f3(( p - f3) - \ oP y(2 - (3 - p)) 

+ 2yp*(l _ /3 )( Ahd _ Aoy ) + 2 yp 2 (l - /3)(A hd - A y) 1 ' 
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We will always assume that the parameters in (|9.13|) lie in the physical region < a, /3,p < 1. 
Under this assumption the two roots of the equation A(y) = lie on the positive real axis, 
so that we may regard $ as defined on two sheets, as we did 0q in the previous section. 



The first sheet corresponds, for y real and small, to J A(y) > in ( |9.13|) . 



From (|9.13|) we see that the singular points of $ (which may coincide for some parameter 
values) are: 

• A simple pole at y = 1. 

• Two square root singularities y±, the roots of the equation A(y) = 0; from (|8.16|) , 

y ± = A± c /A . (9.14) 

Since A < A~ c , these singularities satisfy 1 < y_ < y + . 

• An apparent simple pole at 

Ahd (3{p - 13) ( K v 

m = 17 = x^WY (9 ' 15) 



However, the numerators of the third and fourth terms in ( |9.13| ) may be equal in magnitude 



and opposite in sign when y = y±, cancelling this singularity; from ( |9.15|) and ( |9.13| ) we find 
that this happens when 



p(l -P)-P(2-P-p) = -p(l - P)sJA{y x ). (9.16) 

A little algebra shows that the squares of the two sides in ( |9.16| ) are equal, so that ( |9.16| ) 
holds on the first sheet, and the pole at y\ is absent there, if p(l — (3) — (3(2 — (3 — p) < 0, 
i.e., if (3 > q. 

• An apparent simple pole at y = 0. From (|9.4|) , however, it follows that $ is regular at 



the origin on the first sheet. This can also be seen directly from ( |9.13|) using wA(0) = 1 



We now analyze this generating function in the various regions of the phase plane of the 
system. 

(i) The maximum current region: q < a and q < (3 (region C in figure 2). 

In this region Ao = A~ c , so that from fl9.14|) the square root singularity ?/_ coincides with 
the pole at y — 1; thus A(y) has a factor (1 — y) and from Q9.ll ), 



A(y) = (l-y)(l-p 2 (\ mc ) 2 y). (9.17) 
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Since q < (3, there is no pole at y\ on the first sheet, and thus y — 1 is the singularity 
closest to the origin and controls the asymptotics of the coefficients p n of We write 
Pn = Pn + Pn + Pr? > where is the contribution from the i th term in (|9.13|) (the fourth 
term is regular at y — 1), and calculate the asymptotic form of each pW j n turn. 

The first term in (|9.13|) is /i(y)(l where /i(y) = (2y-l + A~ c py)/ (2y(l+pA~ 
Thus, from flOp and Q , 



/ x (l) + ( s -») = 7: + o(0 



(9.18) 



for any s > 1. The second term is —f2(y)(l — y) 1 ^ 2 , where f2(y) = \Jl — P 2 (^ m c) 2 y / (%(! + 
pA mc )) and we have used Q9.17D . From ([ATP , flOD , and flA"H) , then, 



(2) 



-/a(l)(l-^ ' 



/2a; 



i 



^/T^q~ 1 p 2 + 6g 2 (l-g) 1 



+ 



2 32g 2 v / l — g n^/wn 

The third term is f 3 (y)(l — y) 1 / 2 , where 



+ 0(n- 5 / 2 ) 
— + 0(n~ 5/2 ) 



(9.19) 



/s(y) 



/3(p-/3)v/1-p 2 (A- c ) 2 |/ 
2yp 2 (l -/?)(! +pA- c )(A hd -A- c y) ■ 



(9.20) 



As above 



p (3) 
in 



VT^P(p-P) 1 



+ 0(n- 5/2 ) 



0(n- 5 / 2 ). 



4(/3 — g) 2 n^Jixn 
Adding ( |9.18|) , ( |9.19|) , and ( |9.21|) gives the density p n to order n~ 3 / 2 : 

i v / r^g i 



(9.21) 



irn 



(9.22) 



+ / p 2 + 6g 2 (l - g) + /tyT^p-/?)^ 1 



32g 2 v^^ 



4(/3 — g) 2 I n^/lm 



+ 0{n 



-5/21 



In the p \ limit, with the scaling a = pa, (3 = pf3, this result corresponds to that 
determined in [13|]. The last coefficient ( |9.22|) is singular on the (3 = q boundary of the 
maximum current region and in particular at p = 1; see the comment following (|8.20|) . 
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(ii) The low density region: a < (3 and a < q (region A in figure 2). 



Here Ao = Aid- Since y± = A^ c /Ao > 1, the square root singularities lie strictly to the 
right of the pole at y — 1. Let us again write p n = p^ + pffl + pffl + p^\ with pW the 
contribution from the i th term in ( J9.13 ). Each of the first two terms in (|9.13|) has a simple 
pole at y — 1, so that from ( |A.7| ), 



2y-l + X ld py - JA(y) 



2y(l+p\ 



id) 



o(yZ n ) 



y=l 



a(l — a) 
p — a 2 



+ o(yZ n ), 



(9.23) 



where we have used y A(l) = (p — 2a + a 2 )/ (p(l — a)) (note that this is the positive square 

root). To go further in the asymptotics we must consider separately two subregions of the 
low density region, and their common boundary (see Figure 2). 

Subregion A I: (3 < q. 

In this subregion A id < Ahd < A mc . Since (3 < q, the pole at y\ = Ahd/Aid > 1 lies 
on the first sheet and satisfies 1 < y± < and thus makes the next contribution to the 



asymptotics beyond (|9.23). Thus 



(3 ((p - (3)y/A{y) + (p - g) - \ lA py{2 - (3 - p)) 
2yp 2 (l+pX )(l-(3)X hd 



Ahd 



+ o(s~ 



y=yi 



(l-a)(p-2{3 + {3 2 ) ( a{p-a){l- (3) ' 



n+l 



+ o(s~ 



(9.24) 



for some s > Am/Am- The asymptotics to order o(s ") are obtained by adding (|9.23| ) and 
(PI): 



_ q(l-q) (l- a )(p-2f3 + (3 2 ) ( a{p - q)(l - g) \ w+1 . 

Pn — _ o /_ _9\/i o\ of- o\fi _s r yy~ 



p — cr 



(9.25) 



Further corrections, which arise from the singularity at ?/_, could be calculated; the leading 
order is 0(yZ n /n 3 ^ 2 ). 

Subregion All: (3 > q. 



Now Am < A mc and 1 < the next contribution to the asymptotics beyond ( |9.23| ) 
comes from the square root singularity at present in the second and third terms of 
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( gig ). Now 

A(y) = (l-y/y^(l-yp 2 \ mc \ lA ) 1 
so that, writing p^ 2 )* for the contribution to from this singularity, we have 



(9.26) 



p! 2) * + p! 3) 



1 + pA M 



x 



l-2/P 2 AmAd P(p-P)y/1 -yp 2 \- c \ld 

2y(l-y) + 2yp2(i_ /g ) (Ahd _ Ald2/ ) 



«/=«/- 



2n\/im 



0{yZ n /n b ' 2 ) 



"(p-")(l-p) 1/4 \Am^( A id-^hd) 1 / a(p-a) 



2(p-« 2 )(A mc - Ai d )(A mc - A hd ) n^rm \ (1 - a)(q + 1 - p) 
1 / a(p — a) \ n \ 



n 2 yjn \ (1 — a)(q + 1 — p) 



(9.27) 



From flpD and (|9~27D 

a(l — a) 



p — ar 



a(p-a)(l-p) l ^^X^ c (X ld - A hd ) 1 / a(p-a) 



+ O 



2(p-a 2 )(A mc - Ai d )(A mc - A hd ) n^fwri \ {\ - a)(q + 1 -p) t 
1 / Qi(p-a) \"\ 



(9.28) 



The AI/AII boundary: (3 = q. 



Here ?/i = y_ and the leading correction to p n beyond the constant term ( |9.23|) is an 
0{n~ l l 2 ) contribution from the third term in ( |9.13| ). From ( |9.23| ) and ( |9.26 ), 



Pn 



a(l — a) 
p — a 2 

a(l — a) 
p — a 2 



P(p-P)>Jl -|/P 2 A- c Aid 
2(l+pAid)yp 2 (l-/?)A h d 

a(p — o)(l — p) 1 / 4 



Vi 



+ o 



Vi 



n\/7m 



(9.29) 



y=yi 



1-/3 



(p — a 2 ) 



1 / a(p — a) 
(3(p-(3)^i{f5 2 (l-a) / 



+ O 



nwirn 



The low-density results ( |9.25|) , (|9.28 ), and ( |9.29|) agree with |H| in the p \ limit, 
(iii) The high density region: (3 < a and (3 < q (region B in figure 2). 
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In this region the generating function is obtained by the substitution Ao = Ahd in (|9.13| ). 
After some tedious algebra, this leads to 



P-P 1 



so that the density is constant: 



p- p 2 1-y'' 

P-P 
" p-p 2 ' 



(9.30) 



(9.31) 



10 Exact Expressions for Finite Systems 

In this section we obtain exact and explicit expressions for the current and density profile 
for finite systems and all values of a, P and p. We shall do this by using the algebraic rules 
( |4.3|) and ft4.4|) . This provides a complementary approach to that of sections || and |9| where 
large N properties are calculated directly. 

Our first task is to evaluate z n = {Wi\K n \Vi). We proceed by writing K n as a sum of 
irreducible (with respect to rule (|4.3| )) strings in the following manner 



r=0 g=0 



(10.1) 



It turns out that a n ^ r is given by the expression 
a, 



"n,r 



E 

t=0 L 



n 
r + t 



n — r — 1 
t 



n + 1 
r + t + 1 



n — r — 2 
t-1 



[l-pf. (10.2) 



with the conventions ^ ^ J = 1 and ^ j = for any integer X . The proof of (|10.2|) is 

left to appendix |B[ Here we check a few simple cases. From (|10.2|) and our conventions for 
the binomial coefficients we find 



Q"n,n—1 



a 



n,n—2 



Q"n,n—3 



n-(l-p) , 
-(1-p) 



n 
2 



n 
3 



n 



n + 1 
2 



(1-p) -(1-p) 5 



(10.3) 
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which yield using (|10.1|) 



K = p+(D 1 + E 1 



p + (l+p)(D 1 + E 1 ) + (D 2 1 + E 1 D 1 + E 2 ) , 

2p -p 2 + (2 +p)(D 1 +Et) + (2 + p)(D 2 + E 1 D 1 + E 2 ) 

+(Df + E X D\ + E\D X + E\) , 



;io.4) 



as can be verified by direct calculation. From ( |10.2| ) we determine an exact expression for 
Zj\f by using the action of D\, E\ (|4.4|): 



N 



Z N — X! aN < 



r=0 



(p(l-(3)/(3Y^-(p(l-a)/g) 
(p{l-(3)/(3)-{p{l-a)/a) 



r+l 



(10.5) 



where, without loss of generality, we have taken (Wi|Vi) = 1. Together with ( |4.29|) , ( |10.5| ) 
yields an exact expression for the current. 



To check the limit p — > we use the identity 



E 



N 
X-i 



M 
Y + i 



N + M 
X + Y 



;io.6) 



in order to obtain 



/ In - r - 1 > 
I n — r , 

r(2n — r — 1)! 
n\{n — r)\ 



This agrees with Eq. 39 of [T3| . 



2n — r — 1 
n — r — 1 



Also consider p — 1, then (|10.2j ) becomes a n r = I | and 



z N 



(1 _ - (1 _ a )a-( N+ V 



(10.7) 



(10.8) 



(l-/3)//3-(l- a )/a 

One can check that this recovers the results ( |7.6j ) and ( |7.17|) of section [7[ 

In order to write down the density profile we use an expression derived in appendix O: 



n-l 



D,K n = (1 - p) £ A(r)iT^ + £ a n , rJ D[ +1 , 

r=0 r=0 



(10.9) 
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where a n ^ r is given by (|10.2|) and 

with the convention 1(0) = 1. It can be checked that 

pA{n — 1) = a n fl for n > . 
Inserting ( |10~9D into flOIp yields 

N-i-l 

pZ N -! + (l-p) M r ) Z N-l-r + Zi-1 E a N -i tr (p(l - /?)//?) 



10.10) 



(10.11) 



(rj)jv = 2jv x 



r=0 



JV-i 

E 

r=0 



r+1 



10.12) 



Expression ( |10.12|) together with ( |10.5| ) gives an exact expression for the density profile 
of parallel updating for all system sizes. Through the mappings of section [| it also provides 
exact expressions for the density profiles of ordered and sublattice parallel updating. 



The Case a = (3 = p 

More can be said in the special case of a = (3 = p where many formulae simplify 
considerably. We take advantage of this to simplify the expression for the density profile 
and to write the two-point correlation functions as a sum of one point correlation functions. 



First we note that in this case ( 10.5| ) simplifies as follows: 



EE 

r=0 t=r 



n 
t 



n — r — 1 
t — r 



n + 1 
t + l 



E 

t=0 
n 

E 

t=0 



n 
t 



n + 1 



n + 1 
t + l 



n 
t-1 



n — r — 2 
t-r-1 



1 n+1 
n + 1 { t 



n + 1 
t + l 



(1 - pY = A{n + 1) 



where A(m) is given by (|10.10|) and we have used 

' / M - r x 



E 

r=0 



t — r 



(r + 1) 



M + 2 
t 



{r + l){l-pY 



(10.13) 



10.14) 



From our mapping of the model onto the ASEP with a second class particle, described in 
section |5|, one can check that ( |10.13|) is precisely formula (4.10) in p4| . 
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We also find using the representation (|6.6| )- (|6.8j ) that 



(D 1 +p)K-K(D 1 +p) = K(E 1 +p)-(E 1 +p)K = D X E X -E X D X = {l- v )\V l ){W 1 \ . (10.15) 
Thus the density profile is given by 

fa)* = (r i+1 ) N + ^^-A(i)A(N-i-l) 

(1 ~ P) Eg A(N - n)A{n) + pA(N) 

A(N + 1) + pA(N) ' { } 



One can also use (|10.15|) to relate the two-point correlation function Q4.32| ) to one-point 
correlations: 

(nil-r^N = ( Ti (l - 7j— i)) jv + * A(N-j + i + l) 

Zj N 

+(1-p)^^EW(^- S h, (10-17) 

Z ^ n=0 

= ^" ^(-^-^(iV-i + i + l) 
l=l+1 

E A(iV - I + 1) 2 2 (-p) n (^)/-2-n^-2-n • (10.18) 



«=i+2 n=0 



11 The density in the bulk 

In this section we combine the information derived from generating functions in sections [8] 
and [5] with the the exact calculations of the preceding section to obtain expressions for the 
bulk density of the system, that is, for the large- iV limit of (t^at at constant 9 = i/N. As 
we will see, this bulk density is constant except on the boundary of the low and high density 
regions, and its value is may be guessed by taking the n — > oo limit in the asymptotics of 
section |], that is, in ( |9.22| ), in ( |9.25| ) and ( |9.28| ), and in ( |9.31| ) (where no limit is needed). 



However, it needs to be shown that this limit indeed gives the correct bulk density, and we 
shall do so here. 

The key to the calculation is to study the difference in densities at successive sites. 
Writing z n = z n (a,/3,p), z* = z n (l,fi,p), and z\ = z n (p,p,p), we have from ( |10.12|) , ( |10.5| ), 
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and fll0l3D , 



{Ti)N= j , (11.1) 

ZN + PZN-1 

and hence the density difference Api = (rj)jv — (Tj_i)jv is given by 

A = [P(l ~ 0)/P](zi-iZ*N-i - z i-2 z *N-i+i) ~ (! - g)4r-i-i^-i , 112 v 



We analyze the asymptotic (N — ► oo) behavior of ( |1 1 . 2| ) in various parts of the phase plane; 



properties of the bulk are obtained by the scaling i = 6N, < 9 < 1, and we will also be 
interested in the transition regions i <C N, N — i <C N. The boundary regions in which i or 
N — i remains finite were investigated in section || 



The point a = (3 = p always lies in the maximum current region, so that by (|8.20 ), z\ 



Cnr 3 ' 2 (A~ c )~ n for large n, where here and below C designates some unspecified constant. 
When (a,(3,p) lies in the maximum current region, so does (l,/3,p), so that z n and z* have 
this same asymptotic form and thus for N, i, and N — i large, Api ~ C 'N 3 ' 2 i~ 3 ' 2 (N — i)~ 3 ^ 2 . 
Thus in the bulk Api = 0(N~ 3 / 2 ) and hence the density in the bulk is constant. Moreover, 

N-i 

(TN-i)N ~ Pbulk = ^Pi (ll- 3 ) 

j=8N 

vanishes as i, N — * oo, so that ( |9.22| ) implies that this bulk density has value 1/2. 

When (a, 7) lies in the low density region and (3 < q (i.e., in subregion I), (l,/3,p) lies 
in the high density region, and thus from ( |8.22| ) and reflection symmetry, z n ~ CA^ n and 

z* n ~ CA^d*; since here Ahd < A~ c , Api ~ (^(Aid/Ahd)^ - ^- As above, this implies that the 
bulk density is constant and equal to its value at the left end of the system, which, from 
( |9.31| ) and the reflection symmetry, is a(l — a)/{p — a 2 ). The argument in subregion II is 
similar, with z* ~ Cn- 3 / 2 (X mc )- n and A Pi ~ C(N - _3/2 (A ld /A mc ) (JV - i} ; the bulk density 
is the same. By reflection symmetry the bulk density in the high density phase is constant 
and equal to (p — 13) /{p — f3 2 ). 

The case a = (3 < q requires special attention. Here a slight extension of the arguments 
of section § shows that z n ~ CnX^ 1 . As in the subregion I case above, z* n ~ and the 

z\ term in ( |11.2|) can be neglected; since Ahd = A^, Api ~ CN^ 1 , so that the density profile 



is linear. The values of (rj)jv at the left and right ends of the system are, from ( p.31| ) and 
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the symmetry, a(l — a)/(p — a 2 ) and (p — 0)/{jp — (3 2 ), respectively. This linear profile may, 
as usual, be interpreted clS du superposition of shocks. 



12 Summary 



In this paper, we have presented an exact solution for the steady state of a simple cellular 
automaton describing traffic flow: the ASEP with parallel (synchronous) updating and open 
boundary conditions. The solution is based on recursion relations in the system size for the 
steady-state weights of the configurations, or, equivalently, on formulae for these weights 
as matrix elements of operators satisfying a quartic algebra. By writing these operators as 
rank four tensors, we were also able to express the relevant physical quantities in terms of 
a simpler matrix algebra in which the operators satisfy quadratic relations. 

We used several different methods to extract explicit expressions for observables from 
the matrix algebra. The first applied when p — 1, in which case a two dimensional represen- 
tation of the quadratic algebra exists; this made it possible to obtain analytic expressions, 
in both the finite and the infinite system, for the current and for one and two point cor- 
relation functions. The results confirmed conjectures in |11], p0| . The two point function 
is particularly interesting, because its oscillating behavior directly reflects the particle-hole 
attraction caused by the parallel updating. For a = (3 (still with p = 1) we obtained also 
closed formulae for the fluctuations in the number of particles (cars) in the system. Second, 
for general p we derived exact formulae, in finite systems, for the current and the one point 
function, and for the two point function in the case a = f3 = p; the method was essentially 
an inductive use of the relations of the matrix algebra. The resulting formulae involve rather 
complicated combinatorial expressions in which it is difficult to take the limit of infinite sys- 
tem size. Third, again for general p, we used the analytic properties of generating functions 
to compute asymptotic expressions for the current (and therefore the phase diagram) and 
for density profiles near the boundaries of the system. Finally, we combined the results of 
the last two methods to determine the density in the bulk. 



Our results confirm the phase diagram conjectured in [11]. It is similar to that of 



random sequential updating [13, 14] : there are three phases and, for example, near the right 
boundary we found exponential decay to the bulk density in the low density phase, algebraic 
decay to the bulk density in the maximum current phase, and a constant density profile in 
the high density phase. As p increases, the portion of the phase plane corresponding to the 
maximum current phase shrinks until at p = 1 only the high and low density phases are 
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present. It would be of interest to see if the phase diagram could be predicted by simple 
physical considerations such as those of . 



By considering mappings of the matrix algebra used here to those applicable to other 
updating schemes, we can directly translate all our results (finite size and asymptotics) for 
the current and the profiles to the case where the update of the ASEP is done in discrete 
time but not simultaneously (specifically, with the ordered sequential update and sublattice 
parallel updates). A similar mapping of algebras shows that our results apply to a sys- 
tem of particles on a ring, with one second class particle, in the grand canonical ensemble. 
Since relatively few exact properties of the discrete time updating schemes were previously 



known — essentially only the asymptotic current |L6[ — we obtain new results for these mod- 
els. For example, we verify all conjectured results in Table I of [JTTJ| (these describe bulk 
properties) and derive new formulae for density profiles both for finite systems and asymp- 
totically. The simple translation rules for the current and one-point functions (independent 
of the system size or any other parameters) are surprising, since the two point function of 
the ASEP with parallel updating is very different from that for the other updating schemes. 
It would be interesting to investigate if similar relations are true not only for the ASEP but 
for other models. 
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A Asymptotics from generating functions 

The asymptotic behavior of the coefficients h n of a generating function h(t) = ^nt n 
can frequently be determined, up to order o(s~ n ), from knowledge of the singularities of h 
in a disk \t\ < r with r > s. We analyze below the two cases of this sort. Recall that for 
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any real 7 and complex number t* we have the following Taylor series, 



00 / + \ n 

(l-t/t*r=E Q 7,n - , a 7 , n = Oin^ 1 ) . (A.l) 

n=0 Vt ' 



For application to sections || and ^| we will need the special cases 



a-i, n = 1 , (A.2) 



?m anwirn \n*\/n 



= ■ <aa) 



Case 1: If the only singularities of h in the disc \z\ < r are simple poles at t\, . . . ,t m and 
Cj = lim^^l — t/ti)h(t), then /i — Z)£iCi(l ~~ t/U) is analytic in |t| < r and hence for 
any s < r we have, from ( [A.l] ) and (|A.2p , 



h n = c x t\ n + ■■■ + c m t^ n + o(s- n ) . (A.5) 

Case 2: Suppose that h(t) has simple poles at t%, with Cj defined as above, as well as a 
power singularity at some point to > 0; we assume that \ti\ < to for i = 1, . . . , 7n and that 
h(t) can be written in the form 

h(t)=g(t)(l-t/t y , (A.6) 
where the only singularities of g in the disk are the poles. Then for any k > 0, 



fcn = £ crff + E Mt+^q n + O(n^- fc - 2 to n ) (A.7) 

i=0 j=0 



where 6 3 - = g {j) (t )(-t ) j To verify (|A~7|) write h n = (2m)- 1 § c h(t)t~ n - 1 dt, where the 



contour C has m + 1 components Co, C\, . . . C m . For z > 1, Cj is a small circle, traced 
clockwise, around the point £j, and gives the term Cjt^™ in ( |A.7|) . C follows the circle 



\t\ = s counter-clockwise from just above to just below the positive real axis, then the real 
axis to to + e (for e very small), then a small circle of radius e clockwise around t = to, then 
the real axis tot — s. In evaluating the integral on Co we choose k so that + 7 > (proving 
the result for such a k proves it also for smaller k) and write g(t) = Y,j=obj{l — t/toY + 
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g k+ i(t)(l — t/t Q ) k+1 . The contribution of § Cq bj(l — t/t ) j+J t n 1 dt is precisely b 
There are two terms in the remaining contribution: the integral over the large circle is 
0(s~ n ), while the integral back and forth over the real axis and around the small circle is, 
by our choice of k, a constant multiple of Jf gk+i(t)(l — t/t ) k+1+,y t~ n ~ 1 dt, which is easily 
estimated by the saddle point method to be 0(n _7 ~ fc ~ 2 tg "). 



B Proof of formulae (Jtraj) and (MM 



Proof of flUTM) : 

In this appendix we shall prove 



K n = J2 an, E Er q D\ 

r=0 9=0 



where a n>r is given by the following expression 



E 

t=0 



11 

n — r — t 



n — r — 1 
t 



n + 1 
n — r — t 



(B.l) 



n — r — 2 
t-l 



(l-p)'.(B.2) 



We first require a preliminary result: 

(A + p)E% = (1 - p) n (A + p) + E (1 - p) m Er m (B.3) 

m=0 

which is easy to prove by induction using ( f4.3|) : for n = 1 one has (Di +p)Ei = (1 —p){D\ + 
p) + ; then assuming (|B.3|) and right multiplying by E\ yields 

n+l 

(a+p)^ 1 = (i-p) n [(i-p)(A+p)+^] + E( 1 -p) m ^r m+1 (B- 4 ) 

m=2 

= (1 - p) n+1 [A + p] + E (1 - p) m E?- m+1 , (B.5) 

m=l 

hence ( |B.3j ) is proven by induction. 

Using (|B.3j ) eventually leads to the following recursion for a n ^ r 

n—r 

a n+liT = a n>r -i + E a n,r+m(l - p) m for l<r<n (B.6) 

m=0 
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On+1,0 = E a n,mP(l - VT 

with boundary condition ao,o = 1. To see this left multiply (|B.1D by K 



(B.7) 
(B.8) 



K n+1 = {E 1 + D l + P )Y,a n ,rY, E l" lD l 

r=0 q=0 



n+l 



r-1 



E a n , r _! E E r {~ q D\ + £ a n , r £(1 - p)'-« [D x + p] £>? 

g=0 r=0 q=0 



r=l 



n 7 — 1 r— g— 1 

+ E«n,rE E (1-P) m ^r ? ^? 
r=l g=0 m=0 



(B.9) 



(B.10) 



where we have relabeled the indices r, q in the first term of (|B.10| ) and used (|B.3|) to generate 
the second two terms. Relying on not a little dexterity in relabeling and manipulating sums 
one can develop the second two terms of ( |B.10| ) as follows 



r=0 q=0 



E V [(1 - p) r p + ^i' 4 " 1 ] + E an,r E [(! " P) r+1 "^l + (1 " P) r " V 



r=0 



r=l <j=l 



and 



r=0 



E <W [(1 - P) r p + D[ +1 \ + E a n , r El 1 " P) r-ff tf? 

r=l g=l 



n n— <? 



E <v [(i - pYp + d\ +1 \ + e E v +g (i - P ) r i>? 

r=0 (?=1 r=0 



Ean/E'E^i-p)^-^? 

r=l <j=0 m=0 

n— 1 n r— m— 1 

= E E E an,r(l-p)^r 9 - m ^f 
m=0 r=m+l <j=0 

n— 1 n— m r— 1 

= E EEvh»(i-p) b ^ 

m=0 r=l g=0 
n r— 1 n— r 

= EEEm™(i-p)T^. 

r=l g=0 m=0 



(B.ll) 



(B.12) 
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When the expressions ( |B.12|) , (|B . 1 1|) are inserted back into ( |B.10|) , the second term in the 
square brackets of ( [B.ll|) becomes the q = r component of the first term of QB.10|) , and after 
relabeling indices the third term of (|B.11| ) becomes the q = r component of ( |B.12|) , leading 
to 



n+l r 
r=l g=0 



n r—1 n—r n 

+ E E E «n,r +W (l - p) m E r r q Dt + Y, P(l - PTdn, 
r=l 5=0 m=0 m=0 



From (|B.13[) one can read off (]B.6|)- 



Now assume that be written as 



t=0 



(B.13) 



(B.14) 



Inserting ( |B.14j ) into ( |B.6| ), ( |B.7| ) and ( |B.8| ) respectively yields 



d n +l,r,t — ^n,r-l,t + E d n ,r+m,t-m for 1 < T < iV 



d 



n+l,0,t 



In order to show that 

dm 



m=0 



t t-1 

^ ] dn,m,t—m ^ ] ^n,m, 
m=0 m=0 



t—l—m 



d n ,n,t — Otfi 



b n,r,t 



n \ I n — r — 1 \ ( n + l \ ( n — r — 2 
r + t j ( t j ~ ( r + t+1 )[ t-1 



(B.15) 

(B.16) 
(B.17) 

fB.18) 



satisfies ( |B.15| )-( |B.17| ) we employ two well known identities 



N-M 

E 

M 



N-i 
M-i 



N + 1 
M 



N- 1 \ N-l 
M + M-l 



(B.19) 
(B.20) 



Using (|BTT9|) yields 



^ ] dn,r+m,t—m 
m=0 



n 
r + t 



n — r 
t 



n + l 
r + t + 1 



n — r — 1 
t - 1 



(B.21) 
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Then flB.20| ) becomes 



m=0 



n + 1 
r + 1 



n — r 
t 



n + 2 
r + t + 1 



n — r — 1 
t- 1 



(B.22) 



which is the expression (B.18) required to satisfy (B.15). Similarly with the aid of ( |B.19| ) 



then repeated use of ( |B.20|) one finds 



t t-i 

^ ] d>n,m,t—m ^ ] G^ra,m,t— 1— m 
m=0 m=0 
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n+ 1 
f + 1 



n — 1 
t-1 



n 
t-1 



n 
t- 1 



+ 



n + 1 W n — 1 
t ) \ t-2 



n+1 
t 



n 
t 



n + 2 
t + 1 



n — 1 
t - 1 



(B.23) 



thus satisfying ( B.16|) when d n ^ rj t is given by (|B.18| ). Finally when the conventions 



X 



1 and 



X 



\/X are imposed, ( |B.18[ ) implies 



n,t 



n+1 
n+l + t 



n+1 
-t-1 



-2 
t - 1 



(B.24) 



thus satisfying (|B.17|) . 



Proof of ftWDti 
Here we prove 



n-l 



D 1 K n = (1 - p) £ A{r)K n ' r + On,rD[ 

r=0 r=0 



r+1 



(B.25) 



First we note 



n-l 



[Ei +p] = (l- p) n [Ei +p) + £ (1 - p) m Dr" 



m=0 



(B.26) 



which is proven in a similar fashion to (|B.3|) . 
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To prove ( |B.25|) by induction, one can check the case n — or n — 1, then right multiply 
the rhs of (|[25|) by K, using (|B]2|) to obtain 

n— 1 n 

D,K n+1 = {l-p)Y,A(r)K n+1 - r + Y,an,rD r + 1 [D 1 + E 1 +p\ 

r=0 r=0 

7i — 1 n n 

= (1 - P ) A(r)K n+1 - r + °>n,rD r + 2 + On,r(l " p)^^ + p) 
r=0 r=0 r=0 

n r 

+ EE«n,r(l-Pm +1 ~ m (B.27) 
r=0 m=0 

The third term of ([B.27] ) becomes using ( |B.7| ) 

ta^ll-pf^^+p) = ^—^-a n+x , Q {E x +p) . (B.28) 

r=0 P 



The fourth term of ( |B.27| ) may be developed as follows 



n n—m 

E E a n , r {l-p) m Dl +1 - m = E E a n , r+m (l -p)™^ 1 

r=0 m=0 m=0 r=0 

n n— r „ n 

E E «n, r+ m(l - p) m D[ +1 = + E K+l,r ~ ^-l] (B.29) 

r=0m=0 P r=l 



where ( |B.6| ), ( |B.7| ) have been used to obtain the final equality. Putting ( |B.27| ), ( |B.2£ ) and 
( gggD together yields 

n-1 1—71 ra 

ZW^ 1 = (l-p)E^W^ n+1_r + ^a„ + i,o^ + a n+1 , OJ Di + E a n+i^i +1 + ^i +2 

r=0 P r=l 

n 71+1 

= (l-p)E ^(r)^ n+1 ' r + E ^n+i,,^^ 1 (B.30) 

r=0 r=0 



which agrees with ( ]B.25|) , thereby proving (|10.9 ) by induction 
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Figure Captions 

Fig. 1 The exact truncated correlation function g(i,j) in the case p — 1, for i = 25 versus 
j. The system size is N = 100. The oscillating curve (squares) was obtained for 
a = (5 = 0.9, the other curve (crosses) for a = (3 = 0.1. 

Fig. 2 (taken from [0): Phase diagram for the ASEP with parallel (synchronous) update 
for p — 0.5. C is the maximum current phase, A and B are the low and high density 
phase, respectively. The straight dashed lines are the boundaries between phase A I 
and A II (B I and B II). The curved dashed line is the line given by (|6.1|) and intersects 
the line a = /3ata = /3=l — ^1 — p = q (see section ||). The inserts show typical 
density profiles in the various phases; note that the profile is qualitatively the same in 
region A I (B I) and in the portion of region A II (B II) below the curved dashed line. 
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Figure 1: The exact truncated correlation function g(i,j) in the case p = 1, for i = 25 versus j. 
The system size is N = 100. The oscillating curve (squares) was obtained for a = (5 = 0.9, the 
other curve (crosses) for a = (3 = 0.1. 
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Figure 2: (taken from [T^j): Phase diagram for the A SEP with parallel (synchronous) update for 
p = 0.5. C is the maximum current phase, A and B are the low and high density phase, respectively. 
The straight dashed lines are the boundaries between phase A I and All (B I and B II). The curved 
dashed line is the line given by \6. \ ) and intersects the line a = (3 at a = (3 = 1 — y/\ — p = q (see 
section ^). The inserts show typical density profiles in the various phases; note that the profile is 
qualitatively the same in region A I (B I) and in the portion of region All (BII) below the curved 
dashed line. 
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